#Project: TV polarization
#Author: Josh Ryan and Jeff Lyons
#Working Date: Winter 23
#Task: Use Imai, Kim, Wang PanelMatch analysis--analysis using ps match, Figure 5
#Related files: chamber_level_forstata_v2_112623.dta is loaded


library(foreign)
library(effects)
library(sciplot)
library(sandwich)
library(lmtest)   
library(BMS)
library(rJava)
library(lattice)
library(plyr)
library(lme4)
library(gmodels)
library(MASS)
library(stringr)
library("XLConnect")
library(reshape2)
library(xlsx)
library(dplyr)
library(doBy)
library(Hmisc)
library(stargazer)
library(xtable)
library(dplyr)
library(xtable)
library(foreign)
library(dotwhisker)#https://cran.r-project.org/web/packages/dotwhisker/vignettes/dotwhisker-vignette.html
library(berryFunctions) #for the insert row function
library(ggplot2)
library(gridExtra)
library(wnominate)
library(tidyr)
library(wnominate)
library(oc)
library(haven)
library(PanelMatch)

dat<-read_dta("chamber_level_forstata_v4_102524.dta")
#####Figure 4

#Have to keep making as dataframe because of the way Panelmatch works
dat<-as.data.frame(dat)

#create state-chamber id that is character to use as labels--also because Panelmatch is picky
dat$state_cham2<-paste(dat$state, dat$Chamber,sep = "-")


#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.latebudget<-dat[c("late_budget", "treatment", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "state_cham")]


### Full Analysis with Plots - Updated for PanelMatch 3.1.1

library(ggplot2)
library(gridExtra)

# ========== LATE BUDGET ANALYSIS ==========
dat.latebudget <- as.data.frame(dat.latebudget)
dat.latebudget$state_cham <- as.integer(dat.latebudget$state_cham)
dat.latebudget$year <- as.integer(dat.latebudget$year)

# Create PanelData object for late budget
panel_data_latebudget <- PanelData(
  panel.data = dat.latebudget,
  unit.id = "state_cham",
  time.id = "year",
  treatment = "treatment",
  outcome = "late_budget"
)

results.ps.latebudget <- PanelMatch(
  panel.data = panel_data_latebudget,
  lag = 4,
  refinement.method = "ps.weight",
  match.missing = TRUE,
  covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) + 
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + 
    I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + 
    I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.latebudget <- results.ps.latebudget$att
PE.results.latebudget <- PanelEstimate(sets = results.ps.latebudget, panel.data = panel_data_latebudget)
summary_matrix <- summary(PE.results.latebudget)
coef <- summary_matrix[, 1]
se <- summary_matrix[, 2]
lower <- summary_matrix[, 3]
upper <- summary_matrix[, 4]
time <- c(0, 1, 2, 3, 4)

results.latebudget <- as.data.frame(cbind(coef, se, lower, upper, time))

plot.latebudget <- ggplot(results.latebudget, aes(time, coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks = seq(-.2, .2, .04)) +
  coord_cartesian(ylim = c(-.2, .2)) +
  scale_x_continuous(breaks = seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black")) +
  theme(axis.line.x = element_line(color = "black", size = .5),
        axis.line.y = element_line(color = "black", size = .5)) +
  geom_hline(yintercept = 0, color = "gray") +
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI") +
  ggtitle("Late Budget") + theme(plot.title = element_text(hjust = 0.5))

plot.latebudget
# ========== KURTOSIS ANALYSIS ==========
dat.kurtosis <- dat[c("budget_kurtosis", "treatment", "veto", "mds1", "mds2", "state_ideo", 
                      "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", 
                      "total_intro", "year", "state_cham")]
dat.kurtosis <- as.data.frame(dat.kurtosis)
dat.kurtosis$state_cham <- as.integer(dat.kurtosis$state_cham)
dat.kurtosis$year <- as.integer(dat.kurtosis$year)

panel_data_kurtosis <- PanelData(
  panel.data = dat.kurtosis,
  unit.id = "state_cham",
  time.id = "year",
  treatment = "treatment",
  outcome = "budget_kurtosis"
)

results.ps.kurtosis <- PanelMatch(
  panel.data = panel_data_kurtosis,
  lag = 4,
  refinement.method = "ps.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) + 
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + 
    I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + 
    I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.kurtosis <- results.ps.kurtosis$att
PE.results.kurtosis <- PanelEstimate(sets = results.ps.kurtosis, panel.data = panel_data_kurtosis)
summary_matrix <- summary(PE.results.kurtosis)
coef <- summary_matrix[, 1]
se <- summary_matrix[, 2]
lower <- summary_matrix[, 3]
upper <- summary_matrix[, 4]
time <- c(0, 1, 2, 3, 4)

results.kurtosis <- as.data.frame(cbind(coef, se, lower, upper, time))
results.kurtosis$time<-as.numeric(as.character(results.kurtosis$time))

plot.kurtosis <- ggplot(results.kurtosis, aes(time, coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks = seq(-.2, .2, .04)) +
  coord_cartesian(ylim = c(-.2, .2)) +
  scale_x_continuous(breaks = seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black")) +
  theme(axis.line.x = element_line(color = "black", size = .5),
        axis.line.y = element_line(color = "black", size = .5)) +
  geom_hline(yintercept = 0, color = "gray") +
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI") +
  ggtitle("Budget Kurtosis") + theme(plot.title = element_text(hjust = 0.5))


# ========== PARTY DIFFERENCES ANALYSIS ==========
dat.diffs <- dat[c("diffs", "treatment", "veto", "mds1", "mds2", "state_ideo", 
                   "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", 
                   "total_intro", "year", "state_cham")]
dat.diffs <- as.data.frame(dat.diffs)
dat.diffs$state_cham <- as.integer(dat.diffs$state_cham)
dat.diffs$year <- as.integer(dat.diffs$year)

panel_data_diffs <- PanelData(
  panel.data = dat.diffs,
  unit.id = "state_cham",
  time.id = "year",
  treatment = "treatment",
  outcome = "diffs"
)

results.ps.diffs <- PanelMatch(
  panel.data = panel_data_diffs,
  lag = 4,
  refinement.method = "ps.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) + 
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + 
    I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + 
    I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.diff <- results.ps.diffs$att
PE.results.diffs <- PanelEstimate(sets = results.ps.diffs, panel.data = panel_data_diffs)
summary_matrix <- summary(PE.results.diffs)
coef <- summary_matrix[, 1]
se <- summary_matrix[, 2]
lower <- summary_matrix[, 3]
upper <- summary_matrix[, 4]
time <- c(0, 1, 2, 3, 4)


results.diffs <- as.data.frame(cbind(coef, se, lower, upper, time))
results.diffs$time<-as.numeric(results.diffs$time)


plot.diffs <- ggplot(results.diffs, aes(time, coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks = seq(-.2, .2, .04)) +
  coord_cartesian(ylim = c(-.2, .2)) +
  scale_x_continuous(breaks = seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black")) +
  theme(axis.line.x = element_line(color = "black", size = .5),
        axis.line.y = element_line(color = "black", size = .5)) +
  geom_hline(yintercept = 0, color = "gray") +
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI") +
  ggtitle("Inter-party Polarization") + theme(plot.title = element_text(hjust = 0.5))

# ========== DEM SD ANALYSIS ==========
dat.demsd <- dat[c("dem_sd", "treatment", "veto", "mds1", "mds2", "state_ideo", 
                   "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", 
                   "total_intro", "year", "state_cham")]
dat.demsd <- as.data.frame(dat.demsd)
dat.demsd$state_cham <- as.integer(dat.demsd$state_cham)
dat.demsd$year <- as.integer(dat.demsd$year)

panel_data_demsd <- PanelData(
  panel.data = dat.demsd,
  unit.id = "state_cham",
  time.id = "year",
  treatment = "treatment",
  outcome = "dem_sd"
)

results.ps.demsd <- PanelMatch(
  panel.data = panel_data_demsd,
  lag = 4,
  refinement.method = "ps.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) + 
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + 
    I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + 
    I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.demsd <- results.ps.demsd$att
PE.results.demsd <- PanelEstimate(sets = results.ps.demsd, panel.data = panel_data_demsd)
summary_matrix <- summary(PE.results.demsd)
coef <- summary_matrix[, 1]
se <- summary_matrix[, 2]
lower <- summary_matrix[, 3]
upper <- summary_matrix[, 4]
time <- c(0, 1, 2, 3, 4)

results.demsd <- as.data.frame(cbind(coef, se, lower, upper, time))
results.demsd$time<-as.numeric(results.demsd$time)


plot.demsd <- ggplot(results.demsd, aes(time, coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks = seq(-.2, .2, .04)) +
  coord_cartesian(ylim = c(-.2, .2)) +
  scale_x_continuous(breaks = seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black")) +
  theme(axis.line.x = element_line(color = "black", size = .5),
        axis.line.y = element_line(color = "black", size = .5)) +
  geom_hline(yintercept = 0, color = "gray") +
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI") +
  ggtitle("Democratic Intra-party Polarization") + theme(plot.title = element_text(hjust = 0.5))

# ========== REP SD ANALYSIS ==========
dat.repsd <- dat[c("rep_sd", "treatment", "veto", "mds1", "mds2", "state_ideo", 
                   "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", 
                   "total_intro", "year", "state_cham")]

dat.repsd <- as.data.frame(dat.repsd)
dat.repsd$state_cham <- as.integer(dat.repsd$state_cham)
dat.repsd$year <- as.integer(dat.repsd$year)

panel_data_repsd <- PanelData(
  panel.data = dat.repsd,
  unit.id = "state_cham",
  time.id = "year",
  treatment = "treatment",
  outcome = "rep_sd"
)

results.ps.repsd <- PanelMatch(
  panel.data = panel_data_repsd,
  lag = 4,
  refinement.method = "ps.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) + 
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + 
    I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + 
    I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.repsd <- results.ps.repsd$att
PE.results.repsd <- PanelEstimate(sets = results.ps.repsd, panel.data = panel_data_repsd)
summary_matrix <- summary(PE.results.repsd)
coef <- summary_matrix[, 1]
se <- summary_matrix[, 2]
lower <- summary_matrix[, 3]
upper <- summary_matrix[, 4]
time <- c(0, 1, 2, 3, 4)

results.repsd <- as.data.frame(cbind(coef, se, lower, upper, time))
results.repsd$tscore <- results.repsd$coef / results.repsd$se

results.repsd$time<-as.numeric(results.repsd$time)


plot.repsd <- ggplot(results.repsd, aes(time, coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks = seq(-.2, .2, .04)) +
  coord_cartesian(ylim = c(-.2, .2)) +
  scale_x_continuous(breaks = seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black")) +
  theme(axis.line.x = element_line(color = "black", size = .5),
        axis.line.y = element_line(color = "black", size = .5)) +
  geom_hline(yintercept = 0, color = "gray") +
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI") +
  ggtitle("Republican Intra-party Polarization") + theme(plot.title = element_text(hjust = 0.5))

# ========== LEGISLATIVE PRODUCTIVITY ANALYSIS ==========
dat.prod <- dat[c("billprop", "treatment", "veto", "mds1", "mds2", "state_ideo", 
                  "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", 
                  "total_intro", "year", "state_cham")]
dat.prod <- as.data.frame(dat.prod)
dat.prod$state_cham <- as.integer(dat.prod$state_cham)
dat.prod$year <- as.integer(dat.prod$year)

panel_data_prod <- PanelData(
  panel.data = dat.prod,
  unit.id = "state_cham",
  time.id = "year",
  treatment = "treatment",
  outcome = "billprop"
)

results.ps.prod <- PanelMatch(
  panel.data = panel_data_prod,
  lag = 4,
  refinement.method = "ps.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) + 
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + 
    I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + 
    I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.prod <- results.ps.prod$att
PE.results.prod <- PanelEstimate(sets = results.ps.prod, panel.data = panel_data_prod)
summary_matrix <- summary(PE.results.prod)
coef <- summary_matrix[, 1]
se <- summary_matrix[, 2]
lower <- summary_matrix[, 3]
upper <- summary_matrix[, 4]
time <- c(0, 1, 2, 3, 4)

results.prod <- as.data.frame(cbind(coef, se, lower, upper, time))

plot.prod <- ggplot(results.prod, aes(time, coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks = seq(-.2, .2, .04)) +
  coord_cartesian(ylim = c(-.2, .2)) +
  scale_x_continuous(breaks = seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black")) +
  theme(axis.line.x = element_line(color = "black", size = .5),
        axis.line.y = element_line(color = "black", size = .5)) +
  geom_hline(yintercept = 0, color = "gray") +
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI") +
  ggtitle("Legislative Productivity") + theme(plot.title = element_text(hjust = 0.5))

# ========== COMBINE ALL PLOTS ==========
pdf("treat_effects_v2.pdf", width = 14, height = 8.5)
all.plots <- grid.arrange(plot.latebudget, plot.kurtosis, plot.diffs,
                          plot.demsd, plot.repsd, plot.prod, ncol = 3)
all.plots
dev.off()



# ========== FREQUENCY DISTRIBUTION PLOTS ==========
par(mfrow = c(2, 3))
plot(mset.latebudget, xlim = c(0, 100), ylim = c(0, 25), main = "Late Budget")
plot(mset.kurtosis, xlim = c(0, 100), ylim = c(0, 25), main = "Budget Kurtosis")
plot(mset.diff, xlim = c(0, 100), ylim = c(0, 25), main = "Inter-Party Polarization")
plot(mset.demsd, xlim = c(0, 100), ylim = c(0, 25), main = "Dem. Intra-Party Polarization")
plot(mset.repsd, xlim = c(0, 100), ylim = c(0, 25), main = "GOP Intra-Party Polarization")
plot(mset.prod, xlim = c(0, 100), ylim = c(0, 25), main = "Leg. Productivity")
